#include <iostream>
#include <omp.h>
#include <ctime>
#include <vector>
using namespace std;

#define NUM_THREADS 8
const int SIZE = 100000000;

double GetNum(vector<double> &vec)
{
	double result = 0.0;
	for (int i = 0; i < SIZE; i++)
		result += vec[i] * vec[i];
	return result;
}

double GetPnum(vector<double> &vec, vector<vector<int> > &matrix)
{
	vector<double> tempvec(SIZE);
	for (int i = 0; i < SIZE; i++)
	{
		double temp = 0.0;
		for (int j = 0; j < SIZE; j++)
			temp += vec[j] * matrix[j][i];
		tempvec[i] = temp;
	}
	double result = 0.0;
	for (int i = 0; i < SIZE; i++)
		result += tempvec[i] * vec[i];
	return result;
}

int main(int argc, char *argv[])
{
	clock_t start, finish;
	double totaltime;
	vector<vector<int> > matrix(SIZE);
	for (int i = 0; i < SIZE; i++)
		matrix[i].resize(SIZE);
	srand((unsigned int)time(NULL));
	for (int i = 0; i < SIZE; i++)
	{
		for (int j = 0; j < SIZE; j++)
		{
			if (i == j)
				matrix[i][j] = i +1;
			else
				matrix[i][j] = i + j;
		}
	}
	vector<int> b(SIZE);
	for (int i = 0; i < SIZE; i++)
		b[i] = 100;
	vector<double> x(SIZE);
	vector<double> r(SIZE);
	vector<double> p(SIZE);
	for (int i = 0; i < SIZE; i++)
		x[i] = 1.0;
	for (int i = 0; i < SIZE; i++)
	{
		double temp = 0.0;
		for (int j = 0; j < SIZE; j++)
			temp += matrix[i][j] * x[j];
		r[i] = b[i] - temp;
	}
	for (int i = 0; i < SIZE; i++)
		p[i] = r[i];
	double alpha, beta;
	start = clock();
	while (1)
	{
		alpha = GetNum(r) / GetPnum(p, matrix);

		double betatemp = GetNum(r);
		omp_set_num_threads(NUM_THREADS);
		#pragma omp parallel
		{
			int id = omp_get_thread_num();
			for (int i = id; i < SIZE; i += NUM_THREADS)
				x[i] += alpha * p[i];

			for (int i = id; i < SIZE; i += NUM_THREADS)
			{
				double temp = 0.0;
				for (int j = 0; j < SIZE; j++)
					temp += matrix[i][j] * p[j];
				r[i] -= alpha * temp;
			}
		}

		if (GetNum(r) < 1.0)
			break;

		beta = GetNum(r) / betatemp;
		for (int i = 0; i < SIZE; i++)
			p[i]  = r[i] + beta * p[i];
	}

	finish = clock();
	totaltime = (double)(finish - start) / CLOCKS_PER_SEC;
	cout << "Executive time is " << totaltime << "seconds." << endl;
}